library(MendelianRandomization)
data_dir = "../../dat"
seq_predictions = read.csv(file.path(data_dir, "means_and_uncertainties.csv"))
seq_predictions = subset(seq_predictions, seq_num == 1)
str(seq_predictions)
'data.frame': 3001 obs. of 6 variables:
$ seq_num : int 1 1 1 1 1 1 1 1 1 1 ...
$ mut : Factor w/ 3000 levels "0A","0C","0G",..: 2680 1 2 3 342 343 344 674 675 676 ...
$ X_pred_mean: num -5.63e-05 -1.27e-06 -8.51e-07 -1.21e-06 5.42e-07 ...
$ X_pred_var : num -1.27e-07 -1.61e-07 -1.61e-07 -1.61e-07 -1.61e-07 ...
$ Y_pred_mean: num -9.79e-04 -1.01e-04 -5.30e-05 -3.94e-05 1.14e-04 ...
$ Y_pred_var : num -0.000496 -0.00049 -0.00049 -0.00049 -0.00049 ...
bx <- unlist(seq_predictions["X_pred_mean"])
bxse <- unlist(seq_predictions["X_pred_var"])
by <- unlist(seq_predictions["Y_pred_mean"])
byse <- unlist(seq_predictions["Y_pred_var"])
MRInputObject <- mr_input(bx = bx,
bxse = bxse,
by = by,
byse = byse)
EggerObject <- mr_egger(MRInputObject, robust = TRUE, alpha = .05)
EggerObject
MR-Egger method
(variants uncorrelated, random-effect model)
Number of Variants = 3001
Robust model used.
------------------------------------------------------------------
------------------------------------------------------------------
Residual Standard Error : 11.490
Heterogeneity test statistic = 56644576.7190 on 2999 degrees of freedom, (p-value = 0.0000)
I^2_GX statistic: 100.0%
mr_plot(MRInputObject)
install.packages("devtools")
Installing package into ‘/home/stephenmalina/R/x86_64-pc-linux-gnu-library/3.6’
(as ‘lib’ is unspecified)
trying URL 'https://cloud.r-project.org/src/contrib/devtools_2.2.1.tar.gz'
Content type 'application/x-gzip' length 372273 bytes (363 KB)
==================================================
downloaded 363 KB
* installing *source* package ‘devtools’ ...
** package ‘devtools’ successfully unpacked and MD5 sums checked
** using staged installation
** R
** inst
** byte-compile and prepare package for lazy loading
** help
*** installing help indices
*** copying figures
** building package indices
** installing vignettes
** testing if installed package can be loaded from temporary location
** testing if installed package can be loaded from final location
** testing if installed package keeps a record of temporary installation path
* DONE (devtools)
The downloaded source packages are in
‘/tmp/RtmpCSK8Wq/downloaded_packages’
library(devtools)
Loading required package: usethis
install_github("WSpiller/RadialMR")
Skipping install of 'RadialMR' from a github remote, the SHA1 (18a7a6e9) has not changed since last install.
Use `force = TRUE` to force installation
plot(by/bx, bx/byse, xlim=c(min((by-2*byse)/bx), max((bx+2*byse)/bx)), ylim=c(0, max(bx/byse)))
for (j in 1:length(bx)) {
lines(c((by[j]-1.96*byse[j])/bx[j], (by[j]+1.96*byse[j])/bx[j]),
c(bx[j]/byse[j], bx[j]/byse[j]))
}
abline(h=0, lwd=1); abline(v=0, lwd=1)

library(RadialMR)
r_input <- format_radial(bx, by, bxse, byse)
Missing SNP IDs; Generating placeholders
funnel_radial(egger_radial(r_input, alpha = .05))
Weights not specified: Adopting modified second-order weights
Radial MR-Egger
Estimate Std. Error t value Pr(>|t|)
(Intercept) 17.428069 1.6839374 10.34959 1.088551e-24
Wj 7.145093 0.3373847 21.17788 6.972256e-93
Residual standard error: 91.317 on 2999 degrees of freedom
F-statistic: 448.5 on 1 and 2999 DF, p-value: <2e-16
Q-Statistic for heterogeneity: 25008282 on 2999 DF , p-value: 0
Outliers detected

LS0tCnRpdGxlOiAiTWVuZGVsaWFuIFJhbmRvbWl6YXRpb24gYW5hbHlzaXMgb2YgRGVlcFNFQSBnZW5lcmF0ZWQgZGF0YSIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQpgYGB7cn0KbGlicmFyeShNZW5kZWxpYW5SYW5kb21pemF0aW9uKQpgYGAKCmBgYHtyfQpkYXRhX2RpciA9ICIuLi8uLi9kYXQiCnNlcV9wcmVkaWN0aW9ucyA9IHJlYWQuY3N2KGZpbGUucGF0aChkYXRhX2RpciwgIm1lYW5zX2FuZF91bmNlcnRhaW50aWVzLmNzdiIpKQpzZXFfcHJlZGljdGlvbnMgPSBzdWJzZXQoc2VxX3ByZWRpY3Rpb25zLCBzZXFfbnVtID09IDEpCmBgYAoKYGBge3J9CnN0cihzZXFfcHJlZGljdGlvbnMpCmBgYApgYGB7cn0KYnggPC0gdW5saXN0KHNlcV9wcmVkaWN0aW9uc1siWF9wcmVkX21lYW4iXSkKYnhzZSA8LSB1bmxpc3Qoc2VxX3ByZWRpY3Rpb25zWyJYX3ByZWRfdmFyIl0pCmJ5IDwtIHVubGlzdChzZXFfcHJlZGljdGlvbnNbIllfcHJlZF9tZWFuIl0pCmJ5c2UgPC0gdW5saXN0KHNlcV9wcmVkaWN0aW9uc1siWV9wcmVkX3ZhciJdKQpgYGAKCmBgYHtyfQpNUklucHV0T2JqZWN0IDwtIG1yX2lucHV0KGJ4ID0gYngsIAogICAgICAgICAgICAgICAgICAgICAgICAgIGJ4c2UgPSBieHNlLCAKICAgICAgICAgICAgICAgICAgICAgICAgICBieSA9IGJ5LCAKICAgICAgICAgICAgICAgICAgICAgICAgICBieXNlID0gYnlzZSkKYGBgCgoKYGBge3J9CkVnZ2VyT2JqZWN0IDwtIG1yX2VnZ2VyKE1SSW5wdXRPYmplY3QsIHJvYnVzdCA9IFRSVUUsIGFscGhhID0gLjA1KQpFZ2dlck9iamVjdApgYGAKCgpgYGB7cn0KbXJfcGxvdChNUklucHV0T2JqZWN0KQpgYGAKYGBge3J9Cmluc3RhbGwucGFja2FnZXMoImRldnRvb2xzIikKYGBgCgpgYGB7cn0KbGlicmFyeShkZXZ0b29scykKaW5zdGFsbF9naXRodWIoIldTcGlsbGVyL1JhZGlhbE1SIikKYGBgCmBgYHtyfQpwbG90KGJ5L2J4LCBieC9ieXNlLCB4bGltPWMobWluKChieS0yKmJ5c2UpL2J4KSwgbWF4KChieCsyKmJ5c2UpL2J4KSksIHlsaW09YygwLCBtYXgoYngvYnlzZSkpKQpmb3IgKGogaW4gMTpsZW5ndGgoYngpKSB7CmxpbmVzKGMoKGJ5W2pdLTEuOTYqYnlzZVtqXSkvYnhbal0sIChieVtqXSsxLjk2KmJ5c2Vbal0pL2J4W2pdKSwKYyhieFtqXS9ieXNlW2pdLCBieFtqXS9ieXNlW2pdKSkKfQphYmxpbmUoaD0wLCBsd2Q9MSk7IGFibGluZSh2PTAsIGx3ZD0xKQpgYGAKYGBge3J9CmxpYnJhcnkoUmFkaWFsTVIpCnJfaW5wdXQgPC0gZm9ybWF0X3JhZGlhbChieCwgYnksIGJ4c2UsIGJ5c2UpCmZ1bm5lbF9yYWRpYWwoZWdnZXJfcmFkaWFsKHJfaW5wdXQsIGFscGhhID0gLjA1KSkKYGBgCgo=